function [Z,W,eta,invlambda,zcidx,newK,etime,fval] = mcmc(Z, W, eta, invlambda, X, y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, tnidx, C, ell, poisstrunc, algtype, loopi)
[D,K] = size(Z);
[Nall,tK] = size(W);
Ntrain = numel(invlambda);
M = size(eta, 1);
if K~=tK || K~=size(eta,2) || Nall~=size(X,2) || D~=size(X,1) || Ntrain~=numel(y) || M~=numel(tnidx)-1
    error('Please check data dimension.');
end

if ~any(algtype == 1:2)
    error('cfg:invalidparam', ...
        ['Unknown ''algtype'': %d\n', ...
         '    1 for ''semi-collapsed Gibbs sampling'';\n', ...
         '    2 for ''semi-collapsed EM-style''.'], algtype);
end

etime = 0;
ellp = 0.5*C*ell;
yp = 0.5*C.*y;
ypetaZtX = zeros(1, Ntrain);

% sampling invlambda
tStart = tic;
for m = 1:M
    nind = tnidx(m):tnidx(m+1)-1;
    ypetaZtX(nind) = yp(nind)'.*(eta(m,:)*Z'*X(:,nind));
    tmu = 1./abs(ellp(nind)'-ypetaZtX(nind));
    indinf = isinf(tmu);
    if any(indinf(:))
        tmu(indinf) = max(tmu(~indinf)).^2; % avoid infinite mu
    end
    if algtype == 2
        invlambda(nind) = tmu;
    else
        invlambda(nind) = invnrnd(tmu, 1);
    end
end
tElapsed = toc(tStart);
fprintf('%.4f[%d]: %s (%.2fs) | invlambda\n', C(1), loopi, ...
    num2str([mean(invlambda(:)),std(invlambda(:),1)], ['%.2f',char(177),'%.2f']), tElapsed);
etime = etime + tElapsed;

% sampling W
tStart = tic;
tmp = sqrt(invsigmaxsqr).*Z;
invU = tmp'*tmp;
invU(1:K+1:K*K) = invU(1:K+1:K*K) + invsigmawsqr;
R = choll(invU);
u = R\(R'\((invsigmaxsqr.*Z')*X));
if algtype == 2
    W = u';
else
    W = (u + R\randn(K, Nall))';
end
tElapsed = toc(tStart);
fval = fobj(Z, W, eta, X, y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, tnidx, C, ell);
fprintf('%.4f[%d]: %s (%.2fs) | W\n', C(1), loopi, ...
    num2str([mean(W(:)),std(W(:),1),fval], ['%.2f',char(177),'%.2f',', %.4f']), tElapsed);
etime = etime + tElapsed;

% sampling Z
tStart = tic;
newks = zeros(1,poisstrunc+1);
mZ = sum(Z);
[detsA, suminvA] = scgibbshelper(invsigmawsqr(ones(size(invsigmaxsqr))), invsigmaxsqr, poisstrunc);
for d = randperm(D)
    zrow = Z(d,:);
    mZ = mZ - zrow;
    xi = zeros(Ntrain, K);
    Zdxit = zeros(1, Ntrain); % 1*Ntrain
    for m = 1:M
        nind = tnidx(m):tnidx(m+1)-1;
        xi(nind, :) = kron(yp(nind).*X(d,nind)', eta(m,:));
        Zdxit(nind) = zrow*xi(nind, :)';
    end
    ZdWt = zrow*W';
    ZdV = [ZdWt, Zdxit];
    ypetaZtX = ypetaZtX - Zdxit;
    zrow_sampler(zrow, mZ./(D-mZ), [W;xi], ...
        -[invsigmaxsqr*X(d,:)';1+(ellp-ypetaZtX').*invlambda], ...
        [invsigmaxsqr(ones(Nall,1),1);invlambda], ...
        ZdV, 1, algtype);
    ZdWt = ZdV(1:Nall);
    Zdxit = ZdV(Nall+1:end);
    Z(d,:) = zrow;
    mZ = mZ + zrow;
    ypetaZtX = ypetaZtX + Zdxit;
    
    % semi-collapsed sampling
    tc1 = invsigmaxsqr.*(X(d,:) - ZdWt);
    loglhdk = 0.5*Nall*((1:poisstrunc)*log(invsigmawsqr) - log(detsA)) ...
        + (0.5*sum(tc1.^2)*ones(1,poisstrunc)).*suminvA;
    b = zeros(1, M);
    tc2 = zeros(1, M);
    for m = 1:M
        nind = tnidx(m):tnidx(m+1)-1;
        b(m) = 0.25*(X(d,nind).^2*(invlambda(nind).*C(nind).^2));
        tc2(m) = X(d,nind)*(yp(nind).*(1+(ellp(nind)-ypetaZtX(nind)')'.*invlambda(nind)));
    end
    [detsB, suminvB] = scgibbshelper(invsigmaetasqr(ones(size(b))), b, poisstrunc);
    loglhdk = loglhdk + 0.5*M*((1:poisstrunc)*log(invsigmaetasqr)) ...
        - 0.5*sum(log(detsB)) + 0.5*((tc2.^2)*suminvB);
    logpk = [0, cumsum(log(alphav/D)-log(1:poisstrunc))] + [0, loglhdk];
    if algtype == 2
        [~, newk] = max(logpk);
        newk = newk - 1;
    else % algtype == 1
        newk = logmnrnd(logpk) - 1;
    end
    if newk > 0
        Z(d,K+1:K+newk) = 1;
        mZ(K+1:K+newk) = 1;
        tmu = tc1'.*(suminvA(:,newk)./newk);
        tmu = tmu(:, ones(1,newk));
        tSigma = scgibbshelperinv(detsA, suminvA, newk);
        if algtype == 2
            W(:,K+1:K+newk) = tmu;
        else % algtype == 1
            W(:,K+1:K+newk) = mvnrnd(tmu, tSigma);
        end
        tmu = tc2'.*(suminvB(:,newk)./newk);
        tmu = tmu(:, ones(1,newk));
        tSigma = scgibbshelperinv(detsB, suminvB, newk);
        if algtype == 2
            eta(:,K+1:K+newk) = tmu;
        else
            eta(:,K+1:K+newk) = mvnrnd(tmu, tSigma);
        end
        for m = 1:M
            nind = tnidx(m):tnidx(m+1)-1;
            ypetaZtX(nind) = ypetaZtX(nind) + sum(eta(m,K+1:K+newk))*(yp(nind)'.*X(d,nind));
        end
        K = K + newk;
    end
    newks(newk+1) = newks(newk+1) + 1;
end
newK = sum(newks(2:end).*(1:poisstrunc));
zcidx = sum(Z(:,1:K-newK))==0;
Z(:,zcidx) = [];
W(:,zcidx) = [];
eta(:,zcidx) = [];
K = size(Z,2);
tElapsed = toc(tStart);
fval = fobj(Z, W, eta, X, y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, tnidx, C, ell);
nfeapu = sum(Z,2);
maxnewkid = find(newks, 1, 'last');
fprintf('%.4f[%d]: %d(%d), %s, [%s][%s], %.4f (%.2fs) | Z\n', C, loopi, ...
    K, maxnewkid-1, num2str([mean(nfeapu),std(nfeapu(:),1)], ['%.2f',char(177),'%.2f']), ...
    num2str(newks(2:maxnewkid)./(sum(newks)-newks(1)), '%.2f '), ...
    num2str(newks(1:maxnewkid)./(sum(newks)-cumsum([0,newks(1:maxnewkid-1)])), '%.2f '), ...
    fval, tElapsed);
etime = etime + tElapsed;

% sampling eta
tStart = tic;
for m = 1:M
    nind = tnidx(m):tnidx(m+1)-1;
    tmp = 0.5*(C(nind).*sqrt(invlambda(nind)));
    XtZ = X(:,nind)'*Z;
    tmp = XtZ.*tmp(:,ones(K,1));
    invB = tmp'*tmp;
    invB(1:K+1:K*K) = invB(1:K+1:K*K) + invsigmaetasqr;
    R = choll(invB);
    b = yp(nind).*(1+ellp(nind).*invlambda(nind));
    b = R\(R'\(XtZ'*b));
    if algtype == 2
        eta(m,:) = b';
    else
        eta(m,:) = (b + R\randn(K, 1))';
    end
end
tElapsed = toc(tStart);
fval = fobj(Z, W, eta, X, y, alphav, invsigmawsqr, invsigmaetasqr, invsigmaxsqr, tnidx, C, ell);
fprintf('%.4f[%d]: %s (%.2fs) | eta\n', C(1), loopi, ...
    num2str([mean(eta(:)),std(eta(:),1),fval], ['%.2f',char(177),'%.2f',', %.4f']), tElapsed);
etime = etime + tElapsed;

end